16. 混合ガウシアンモデル(Gaussian Mixture Model, GMM)#
16.1. EMアルゴリズム#
EMアルゴリズムは、EステップとMステップを交互に繰り返すことで、隠れ変数を含む確率モデルのパラメータを推定する手法である。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# ----------------------------
# データ生成 (Ground Truth)
# ----------------------------
np.random.seed(42)
true_means = np.array([[0, 0], [3, 3], [-3, 3]])
true_covs = [
np.array([[1, 0.3], [0.3, 1]]),
np.array([[1, -0.2], [-0.2, 1.5]]),
np.array([[1.5, 0.4], [0.4, 1]])
]
true_weights = np.array([0.3, 0.4, 0.3])
n_samples = 10**3
X = np.vstack([
np.random.multivariate_normal(mean, cov, int(w * n_samples))
for mean, cov, w in zip(true_means, true_covs, true_weights)
])
# ----------------------------
# ガウス分布のPDF
# ----------------------------
def gaussian_pdf(x, mean, cov):
d = x.shape[1]
det = np.linalg.det(cov)
inv = np.linalg.inv(cov)
norm_const = 1.0 / np.sqrt((2 * np.pi) ** d * det)
diff = x - mean
return norm_const * np.exp(-0.5 * np.sum(diff @ inv * diff, axis=1))
# ----------------------------
# プロット: データと推定されたガウス分布の等高線
# ----------------------------
def plot_gaussian_ellipse(mean, cov, ax, color='red'):
vals, vecs = np.linalg.eigh(cov)
order = vals.argsort()[::-1]
vals, vecs = vals[order], vecs[:, order]
angle = np.degrees(np.arctan2(*vecs[:, 0][::-1]))
width, height = 2 * np.sqrt(vals)
ellip = Ellipse(xy=mean, width=width, height=height,
angle=angle, edgecolor=color, fc='None', lw=2)
ax.add_artist(ellip)
def snapshot_gaussians(X, mu, sigma, K, iteration):
plt.figure(figsize=(5, 5))
plt.scatter(X[:, 0], X[:, 1], s=10, alpha=0.5)
colors = ["red", "green", "blue"]
for j in range(K):
plot_gaussian_ellipse(mu[j], sigma[j], plt.gca(), color=colors[j])
plt.scatter(mu[j, 0], mu[j, 1], c=colors[j], marker="x", s=100)
plt.title(f"Data and Fitted Gaussian Components @ Iteration = {iteration}")
plt.xlabel("x1")
plt.ylabel("x2")
plt.show()
# ----------------------------
# EMアルゴリズム (GMM)
# ----------------------------
def EM_GMM(X, K, num_iters=100, tol=1e-5):
n, d = X.shape
np.random.seed(11)
# 初期化
phi = np.ones(K) / K
mu = X[np.random.choice(n, K, replace=False)]
sigma = np.array([np.cov(X, rowvar=False) for _ in range(K)])
log_likelihoods = []
snapshots = []
for iteration in range(num_iters):
# E-step
responsibilities = np.zeros((n, K))
for j in range(K):
responsibilities[:, j] = phi[j] * gaussian_pdf(X, mu[j], sigma[j])
responsibilities /= responsibilities.sum(axis=1, keepdims=True)
# M-step
N_k = responsibilities.sum(axis=0)
phi = N_k / n
mu = (responsibilities.T @ X) / N_k[:, np.newaxis]
for j in range(K):
diff = X - mu[j]
sigma[j] = (responsibilities[:, j][:, np.newaxis] * diff).T @ diff / N_k[j]
# 対数尤度
ll = np.sum(np.log(np.sum([
phi[j] * gaussian_pdf(X, mu[j], sigma[j]) for j in range(K)
], axis=0)))
log_likelihoods.append(ll)
# スナップショット保存
snapshots.append((mu.copy(), sigma.copy(), iteration))
# 収束判定
if iteration > 0 and abs(log_likelihoods[-1] - log_likelihoods[-2]) < tol:
print(f"Converged at iteration {iteration}")
break
# アニメーション作成
fig, ax = plt.subplots(figsize=(5, 5))
scatter = ax.scatter(X[:, 0], X[:, 1], s=10, alpha=0.5)
title = ax.set_title("")
def animate(i):
ax.clear()
ax.scatter(X[:, 0], X[:, 1], s=10, alpha=0.5)
mu_snap, sigma_snap, iter_snap = snapshots[i]
colors = ["red", "green", "blue"]
for j in range(K):
plot_gaussian_ellipse(mu_snap[j], sigma_snap[j], ax, color=colors[j])
ax.scatter(mu_snap[j, 0], mu_snap[j, 1], c=colors[j], marker="x", s=100)
ax.set_title(f"Data and Fitted Gaussian Components @ Iteration = {iter_snap}")
ax.set_xlabel("x1")
ax.set_ylabel("x2")
anim = FuncAnimation(fig, animate, frames=len(snapshots), interval=500)
plt.close(fig)
display(HTML(anim.to_jshtml()))
return phi, mu, sigma, responsibilities, log_likelihoods
# ----------------------------
# 実行
# ----------------------------
K = 3
phi, mu, sigma, responsibilities, log_likelihoods = EM_GMM(X, K)
print("Estimated mixing coefficients (phi):", phi)
print("Estimated means (mu):\n", mu)
print("Estimated covariances (sigma):\n", sigma)
# ----------------------------
# プロット: 対数尤度の推移
# ----------------------------
plt.figure(figsize=(6, 4))
plt.plot(log_likelihoods, marker="o")
plt.title("Log-Likelihood over EM Iterations")
plt.xlabel("Iteration")
plt.ylabel("Log-Likelihood")
plt.grid(True)
plt.show()
Converged at iteration 57
Estimated mixing coefficients (phi): [0.3103679 0.30013538 0.38949671]
Estimated means (mu):
[[ 0.08871446 0.04564594]
[-3.08442868 3.07309071]
[ 3.05915799 3.15711097]]
Estimated covariances (sigma):
[[[ 1.01691342 0.38486941]
[ 0.38486941 0.98651172]]
[[ 1.44451192 0.28724678]
[ 0.28724678 0.91035083]]
[[ 1.09349783 -0.24274844]
[-0.24274844 1.27556857]]]
true_means = np.array([[0, 0], [3, 3], [-3, 3]])
true_covs = [
np.array([[1, 0.3], [0.3, 1]]),
np.array([[1, -0.2], [-0.2, 1.5]]),
np.array([[1.5, 0.4], [0.4, 1]])
]
true_weights = np.array([0.3, 0.4, 0.3])